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SECTION  I 
INTRODUCTION 

In  recent  years,  increasing  use  has  been  made  of  the  Doubly  Asymptotic  Approximation 
(DAA)  [1-3)  for  the  transient  response  analysis  of  submerged  structures.  Application  of  the 
DAA  in  conjunction  with  discrete  methods  of  structural  analysis  leads  to  a large  set  of 
coupled  ordinary  differential  equations  (ODE)  in  time.  Such  equations  are  either  solved 
directly  in  terms  of  the  physical  coordinates,  or  indirectly  through  a dimensionality- 
reduction  transformation  from  physical  to  generalized  coordinates.  This  report  deals  exclu- 
sively with  the  first  approach. 

The  physical  coordinates  associated  with  the  direct  time  integration  approach  are  the 
response  degrees-of-f reedom  of  the  structural  model  (the  structural  displacements)  and  the 
nodal  values  of  the  scattered  pressure  field  at  the  contact  surface  (the  fluid  pressures). 
Fluid-structure  interaction  effects  result  in  the  coupling  of  the  fluid  pressures  and  struc- 
tural displacements  normal  to  the  contact  surface  (the  "wet  displacements")  through  a matrix 
differential  equation. 

1.1  SOLUTION  PROCEDURES 

Three  general  schemes  can  be  followed  for  organizing  the  response  calculations: 

• Pressure  Elimination  Solution.  Elimination  of  the  fluid  pressure  equations  results 
in  a differential  system  of  equations  in  the  structural  displacements  only;  this 
system  is  of  third  order  in  time. 

• Simultaneous  Solution.  The  fully  coupled  system  is  processed  as  one  enticy,  i.e., 
both  displacements  and  pressures  are  simultaneously  advanced  through  the  use  of  a 
fully-impllcit  integration  scheme. 

• Staggered  Solution.  The  integration  process  is  carried  out  in  alternating  stages: 
the  solution  of  either  the  structural  or  fluid  system  is  advanced  first,  then  the 
other  system  is  solved  following  a time  extrapolation  of  the  coupling  terra(s)  as  a 
forcing  function. 

The  first  approach  has  the  advantage  of  reducing  the  number  of  equations  to  be  solved, 
but  introduces  unsymmetric  coefficient  matrices  that  are  fully  populated  in  entry  positions 
pertaining  to  the  "wet"  degrees  of  freedom.  Furthermore,  the  appearance  of  third-order  (or 
higher)  temporal  derivatives  of  the  displacement  coordinates  causes  numerical  difficulties  in 
the  treatment  of  initial  and  jump  conditions  in  shock-excited  problems.  (These  difficulties 
can  be  eliminated  by  passing  to  integro-dif ferential  forms,  but  then  time  integrals  of 
applied  forcing  terms  must  be  carried  along  in  the  calculations.) 

The  second  approach  can  be  organized  in  terms  of  symmetric  matrices  and  does  not  re- 
quire special  treatment  to  account  for  acceleration-derivative  terms.  Processing  times 
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ror  large-scale  problems  become  rapidly  prohibitive,  even  in  iinear  analysis,  on  account  of 
the  matrix  connectivity  introduced  by  pressure/wet-displacement  coupling  blocks  that  may 
extend  across  thousands  of  degrees  of  freedom.  For  example,  the  initial  factorization  of  the 
dynamic  coefficient  matrix  for  a model  with  5000  structural  and  250  fluid  equations  (a  typ- 
ical mix)  would  require  roughly  3 hrs  of  CPU  time  on  the  CDC  6600  computer. 

1.2  STAGGERED  SOLUTION  PROCEDURE 

The  staggered  solution  scheme  consists  of  splitting  the  time  integration  task  between 
two  loosely  coupled  processors,  a fluid  and  a structural  analyzer,  with  the  interaction 
effects  being  incorporated  through  an  extrapolation  mechanism.  If  certain  practical  restric- 
tions (described  in  the  body  of  the  paper)  are  met,  this  strategy  offers  the  important  ad- 
vantage of  preserved  program  modularity , in  the  sense  that  the  organization  of  the  structural 
time  integration  package  is  not  sensibly  affected.  A general-purpose  fluid  analyzer  may  be 
developed  and  checked  out  as  independent  of  any  specific  structural  code  and  eventually 
"plugged  in"  as  a modular  component  to  interface  with  any  large-scale  structural  analyzer. 

Furthermore,  subsequent  improvements  made  in  the  fluid  analyze"  will  be  essentially  trans- 
parent to  the  user,  resulting  in  corresponding  gains  in  efficiency.  Additional  advantages 
accruing  in  the  treatment  of  nonlinear  structural  response  problems  are  discussed  in  the 
section  on  implementation. 

Given  the  usual  proportion  of  structural  to  fluid  equations  (4:1  through  25:1  in  three- 
dimensional  problems),  the  bulk  of  the  computational  effort  is  most  likely  to  fall  upon  the 
structural  analyzer,  which  perceives  the  fluid  only  as  an  external  force  environment.  It  is 
therefore  reasonable  to  expect  that  the  fluid-structure  problem  can  be  processed  for  only  a 
marginal  increment  in  the  cost  required  for  processing  the  "dry"  structure  only  (provided, 
of  course,  that  similar  time  increments  and  structural  solution  strategies  are  used  in  both 
cases).  The  price  paid  for  these  computational  advantages  is  the  fact  that  satisfactory 
numerical  stability  characteristics  are  much  harder  to  achieve  for  the  staggered  procedure. 

No  such  difficulties  arise  in  the  pressure-elimination  or  simultaneous  solution  procedures, 
for  which  the  selection  of  one  of  the  many  available  A-stable  time  integration  operators 
suffices  to  secure  unconditional  stability. 

The  purpose  of  this  report  is  to  examine  the  fundamental  algorithmic  properties  — 
convergence,  stability  and  accuracy  — of  the  staggered  solution  procedure  when  the  fluid- 
structure  interaction  is  treated  by  the  DAA.  It  will  be  shown  that  the  simplest  (conven- 
tional) formulation  of  the  staggered  procedure  for  the  structure/DAA  equations  suffers  from 
severe  time-increment  limitations  due  to  stability  considerations.  Such  restrictions 
effectively  rule  out  the  use  of  the  conventional  staggered  solution  procedure  as  a general- 
puroose  method,  although  that  scheme  can  work  satisfactorily  in  a restricted  class  of 
shock-excited  problems. 
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The  main  objective  of  this  investigation  is  to  find  ways  to  extend  the  range  of 
applicability  of  the  staggered  solution  procedure.  It  happens  that  a satisfactory  stabiliza- 
tion of  that  strategy  can  be  achieved  by  a judicious  modification  of  the  original  equations 
of  motion  through  a process  of  augmentat ion . There  are  a number  of  ways  in  which  the  fluid 
and/or  the  structural  equations  can  be  augmented  to  produce  the  desired  stabilization,  but 
only  one  survives  if  the  practical  constraint  of  limiting  modifications  to  existing  struc- 
tural analyzers  is  imposed. 

1.3  OUTLINE  OF  THIS  REPORT 

In  the  section  following,  equations  of  motion  that  govern  a spatially-discretized, 
fluid-structure  interaction  problem  treated  by  the  DAA  are  introduced.  An  associated  two 
degree-of-f reedom  system  is  derived  and  converted  to  dimensionless  coordinates.  This  system 
is  time-discretized  through  application  of  linear  multistep  operators  and  the  resulting 
integration  process  is  expressed  in  the  conventional  staggered  solution  form.  A study  of  the 
domain  of  convergence  of  the  iterated  process  and  of  the  temporal  stability  characteristics 
shows  that  this  formulation  has  a limited  applicability  range  because  of  time-increment 
limitations. 

The  source  of  instability  of  the  conventional  staggered  procedure  is  then  explored  in 
detail.  This  investigation  is  carried  out  by  using  a time-delayed  continuous  form  of  the 
extrapolated  coupling  term  (predictor)  for  deriving  a difference-differential  system.  Exam- 
ination of  the  characteristic  equation  of  this  system  clearly  shows  that  the  onset  of  in- 
stability is  caused  by  the  delayed  feedback  of  fluid  radiation  damping  from  the  fluid 
equation  into  the  structural  equation.  This  result  leads  us  to  consider  the  use  of  "damping 
augmentation"  strategies  that  are  often  employed  in  the  field  of  control  theory  to  stabilize 
systems  with  dead-time  (delay)  elements. 

Several  stabilized  formulations  are  then  developed  by  tailoring  the  governing  equations 
of  motion  in  such  a way  that  appropriate  damping  terms  appear  in  the  structural  equations, 
the  fluid  equations,  or  both.  A study  of  the  iterative  convergence  and  temporal  stability 
of  these  forms  show  that  they  are  globally  convergent  and  remain  unconditionally  stable  for 
certain  extrapolators . 

An  extensive  series  of  numerical  experiments  has  been  conducted  with  the  dual  objective 
of  verifying  the  predictions  of  the  stability  analysis  and  of  assessing  the  global  accuracy 
of  the  computed  solutions  oi  the  two  degree-of-f reedom  system.  Representative  results  of 
this  series  are  presented.  Practical  considerations  regarding  the  computer  implementation 
of  the  stabilized  formulations  to  treat  large-scale  problems  are  then  offered;  special  em- 
phasis is  placed  on  software  modularity  requirements.  Finally,  the  main  conclusions  derived 
from  the  present  investigation  are  summarized. 


SECTION  II 


GOVERNING  EQUATIONS 

Consider  a structure  interacting  with  an  acoustic  medium  through  a contact  boundary  B, 
henceforth  referred  to  as  the  "wet  surface."  The  structure  and  fluid  are  spatially  dis- 
cretized through  the  application  of  finite-element  and  boundary-integral  techniques,  respec- 
tively. (It  is  important  to  note  that  the  corresponding  meshes  on  B are  not  necessarily 
identical.)  The  resulting  matrix  equations  of  motion  may  be  written  as 

Mil+Du  + Ku  = fs  + r - T A (ES  + j,1) 

(1) 

M f £ + p c A £ = p C M f (TT  u - u1) 

where  the  first  set  of  equations  expresses  dynamic  equilibrium  in  terms  of  the  structural 
displacements,  and  the  second  set,  the  DAA,  describes  the  f xuid-structure  interaction  at  the 
wet  surface.  (Recall  that  the  number  of  structural  equations  is  usually  much  greater  than 
the  number  of  fluid  equations.)  In  (1),  M,  I),  and  K are  the  structural  ma  ■ , damping 
and  linear  (or  linearized)  stiffness  matrices,  respectively;  £,  f^  and  £ are  the  struc- 
ture response  displacements,  dry-structure  applied  force,  and  nonlinear  residual  (pseudo- 

S • I 

force)  vectors,  respectively;  £ = £ is  the  scattered  pressure  vector  and  £ is  the 

incident  pressure  vector,  appropriate  to  the  fluid  grid  on  B ; A is  a diagonal  matrix 
embodying  elemental  areas  of  the  fluid  mesh  on  the  wet  surface  B,  Mj.  is  the  fluid  mass 
matrix  as  determined  from  an  analysis  of  incompressible  fluid  motion  appropriate  to  a dis- 
tribution of  elemental  sources  on  B,  T is  a generally  rectangular  transformation  matrix 
that  relates  structural  displacements  to  the  control  points  of  the  fluid  grid  on  B , p 
and  c are  fluid  density  and  speed  of  sound,  respectively,  and  superscript  T denotes 
matrix  transposition;  finally,  the  superscript  dot  (')  denotes  temporal  differentiation. 

2.1  THE  MODEL  PROBLEM 

The  homogeneous  linear  (or  linearized)  part  of  (1)  can  be  expressed  as 


I. 


! 


f 


I 

{ 

I 
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In  the  following  study,  the  structural  damping  term  Du  will  be  systematically  neglected, 
as  its  effect  on  the  response  is  in  most  cases  negligible  when  compared  to  that  of  the  fluid 
radiation  damping  term  TA j.  An  appropriate  two  degree-of-f reedom  problem  associated 
with  the  system  (2),  after  setting  D = 0 , is 


m w + k w = - a g 
s s 

(3) 

mf  g + pcag  = pcmfw 


where  mg,  kg,  a and  m^.  are  generalized  quantities  resulting  from  the  projection  of  M, 
K,  A and  M^.,  respectively,  on  normal  coordinates  w and  g that  diagonalize  the  left- 
hand  side  of  (2).  The  derivation  of  (3)  is  presented  in  Appendix  A.  We  note  that  (3a)  rep- 
resents a pressure-excited  undamped  mechanical  oscillator  whereas  (3b)  represents  a velocity- 
excited  pressure-decay  equation. 

The  system  (3)  can  be  further  reduced  to  the  non-dimensional  form 


£ x + to 


- y 


y + p y = 


(A) 


through  introduction  of  the  dimensionless  variables 


x = w/Z  , y = g/(pcZ) 


2 2 2 

£ = m / (pZa)  , ii  = k l /(£m  c ) 
s s s 


p = mg/(£mf)  = p!la/mf 


(5) 


( ‘ ) = 3/3t  , t = c t/Z 


la  Eqs.  (5),  J. 
merged  cylinder 
fluid  mass),  in 
exponent.  Note 
with  respect  to 


denotes  a characteristic  length  of  the  problem,  e.g.,  the  radius  of  a sub- 
or  sphere;  £ is  a "buoyancy  ratio"  (structural  mass  divided  by  displaced 
is  a reduced  vibration  frequency,  and  p is  a generalized-pressure  decay 
that  the  dot  superscript  has  been  redefined  in  (5d)  to  denote  differentiation 
the  dimensionless  time  t , rather  than  physical  time  t. 
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2.2  TEMPORAL  DISCRETIZATION 

We  shall  consider  the  use  of  implicit,  one-derivative  linear  multistep  (LMS)  methods 
to  effect  the  time  discretization  of  (4).  For  a constant  step-size  h = Ax  , an  m-step 
method  can  be  presented  in  the  compact  form 

z = <5  i + hZ  (6) 

n n n 

where  z stands  for  a scalar  or  vector  state  variable,  the  subscript  n is  the  time  station 

index,  6 is  the  generalized  time  step  8h  (8  being  a method  coefficient),  and  the  his- 
2 

torical  vector  is  a linear  combination  of  m past  solutions: 

m 

hZ  = -T,  (a.  z . - 6.  h i .)  (7) 

n f-*.  l n-i  l n-i 

1=1 

Specific  LMS  integrators  are  characterized  by  the  coefficients  8 , and  6^.  For  example, 
the  trapezoidal  rule  (m  = 1)  has  = -1  , 8 = 6j  = 0.5. 

For  reasons  to  be  justified  later,  we  introduce  a pair  of  LMS  integrators 

x = 6 x + hX  , 6 =6h 

n x n n x x 

(8) 

y =6y  + hy  , 6 =8h 
n y n n y y 

to  treat  the  reduced  equations  (4).  The  resulting  algebraic  system  is 

a + <52w2)  x = -«2y  + C (hx  + 6 hX) 

x n x n n x n 

(9) 

(1+6  p)y  =x  -phy 
y n n n 

The  secondary  state  variables,  velocity  x^  and  pressure  integral  y^  , are  calculated  from 


the  integrators  (8),  i.e. 


x =(x  -h)/6 

n n n x 


y = h^  + 6 y 
n n y n 


Eqs.  (9)  and  (10)  form  the  basis  for  the  analysis  of  the  conventional  staggered  solution 
procedure. 


2.3  RANGE  OF  NON-DIMENSIONAL  PARAMETERS 

The  time-discretized  system  (9)  contains  four  non-dimensional  parameters:  £ , m , p 
and  h . The  first  three  embody  spatial  characteristics  of  the  problem  whereas  h intro- 
duces the  effect  of  the  time  integration  stepsize.  It  is  of  interest  to  the  forthcoming 
discussions  on  the  applicability  of  various  implementations  of  the  staggered  procedure,  to 
exhibit  typical  ranges  assumed  by  such  quantities.  This  information  is  collected  in  Table  1. 

For  £ , u)  , and  p , two  limit  conditions  and  an  illustrative  case  are  shown.  The 

cavity  condition  is  the  limit  of  modal  motions  heavily  dominated  by  the  inertia  of  the  fluid, 

viz.  low-frequency  motions  of  a very  thin  submerged  shell.  The  dry  mode  condition  is 
realized  by  structural  motions  that  do  not  interact  with  the  fluid,  such  as  torsional  modes 
of  structures  of  revolution.  The  illustrative  problem  of  a submerged  spherical  shell  is  of 

interest  because  this  is  one  of  the  few  geometries  amenable  to  exact  analytical  modal  treat- 

ment, as  discussed  in  Appendix  B.  Ranges  quoted  for  the  dimensionless  time  increment  h = 
cAt/S.  are  tabulated  according  to  the  temporal  characteristics  of  the  excitation,  which 
determines  the  energy-spectrum  characteristics  of  the  response. 
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TABLE  1.  Range  of  Non-Dimensional  Parameters 


Parameter 

Limit 

Cases 

Illustrative  Case 

Cavity 

\-°*  ks~° 

Dry  Structural  Mode 

mj.~*0,  a— 0 

Submerged  Spherical 

Shell3 

U) 

0b 

20 

L 2 

0 to  ~ pee  n 

£ 

0 

oo 

2 

~ £e  to  ~ pen 

V 

1 

indet . 

1 to  n + 1 

Q 

Shock-Excited  Problems 

Structural  Dynamics 

Early-Time 

Late-Time° 

Problems  ^ 

Response 

Response 

h 

0.01-0. I 

1-100 

»1 

=Ps/p  = ratio  of  shell  and  fluid  densities,  c = c /c  = ratio  of  sound  speeds  in 
shell  and  fluid,  e = thickness-to-radius  ratio,  n = highest  circumferential  wave 
number  retained.  These  equations  are  derived  in  Appendix  B,  where  numerical  values 
of  to  , £ and  P are  tabulated  for  a particular  combination  of  p , c and  e. 

b k tends  to  zero  faster  than  m . 
s s 

c Problems  characterized  by  wave  propagation  effects. 

^ Period  during  which  shock  wave  starts  to  envelop  the  structure  (t  < 1); 
characterized  by  high-frequency  structural  motions,  high  radiation  damping, 
relatively  small  hydrodynamic  inertia  forces. 


Period  characterized  by  low-frequency  structural  motions,  dominant  hydrodynamic 
inertia,  and  low  radiation  damping  (usually  t ? 10). 


f 


Problems  characterized  by  low-frequency  motions  throughout. 


SECTION  III 


PRESSURE  EXTRAPOLATION  FORMULATION 


The  simplest  formulation  of  the  staggered  scheme  suggested  by  the  format  of  Eqs.  (9),  is 
the  pressure  extrapolation  (PE)  formulation,  which  may  be  described  as  follows.  Assume  that 

. p 

solutions  up  to  the  (n-l)th  time  station  are  known.  We  first  predict  the  pressure  y^  at 
t^  , insert  this  into  the  right-hand  side  of (9a),  solve  for  the  structural  displacement  x^  , 
obtain  the  velocity  from(lOa),  insert  this  into  (10b)  and  finally  solve  for  the  pressure 

yn  . This  corrected  value  may  be  used  in  an  iterative  setting  if  desired.  We  proceed  now  to 
analyze  the  iteration  convergence  and  temporal  stability  properties  of  the  PE  formulation. 

3.1  ITERATION  CONVERGENCE 

Using  k as  an  iteration  cycle  index,  the  iterated  PE  scheme  may  be  written 

E x(k)  = -6  hx  y(k_1)  + b 
x n x n x 

k = 1,2  ...  (11) 

E y(k)  = (x(k)  - hX)/6  + b 
y n n n x y 


2 2 2 2 2 
e - i + e n , n - « h /e 

X X 


1 + ey  T , T = p h 


X = h/£ 


b = hX  + 6 hX  , b = -p  hy 
x n x n y n 


The  initial  y^^  = y^  may  be  obtained  from  a simple  predictor  such  as 

yn0)  = (1  + Y)  Vl  - Y yn-2  (13) 

where  y is  an  extrapolation  parameter  (actually,  extrapolation  corresponds  to  y > 0 . and 
interpolation  to  y « 0)  . 

The  iterative  scheme  (11)  can  be  recast  into  the  standard  form 

z(k)  = R z(k_1)  + b (14) 
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where  z = (x  , y ) , b embodies  terms  independent  of  k , and  the  iteration  matrix  R 
— n n — ~ 


R = 


E 0 
x 


-1 


o e hx 

x 


0 0 


*x  hx 

E E 
x y 


0 E 


0 1 


(15) 


The  iteration  (14)  converges  if  and  only  if  the  spectral  radius  k of  R (largest  eigen- 
value modulus)  is  less  than  1.  Since  k(R)  = 6 \/ (E  E ) , the  convergence  condition  is 

***  x x y 


h/c  < (i  + ex  a ) (i  + sy  'R) /ex 


(16) 


The  right-hand  side  of  (16)  is  minimized  for  modal  motions  with  ft  -*•  0 , f -*•  0 , in  which 
reads 


X < l/6x  (17) 

The  structural  equations  will  normally  be  treated  with  an  A-stable  integrator  to  account  for 
the  wide  frequency  spectrum  present  in  most  discrete  models.  But  for  all  A-stable  LMS 
methods,  8x  lies  in  the  range  of  ^ to  1.  Consequently, 

X = h/5  < 2.0  (18) 

is  the  best  that  can  be  achieved  by  using  the  trapezoidal  rule  (8  = h)  in  (8) . Note  that 

the  convergence  conditions  are  not  only  independent  of  the  pressure  predictor  (as  one  would 
expect  of  a linear  system),  but  of  the  historic  term  composition  of  the  integrator  pair  (8) 
as  well, 

3.2  TEMPORAL  STABILITY 

The  numerical  stability  of  the  PE  formulation  has  been  studied  for  a variety  of  LMS 
integrators,  assuming  the  two-step  predictor  (13),  and  a fixed  number  (k)  of  passes  per 
time  step.  No  details  of  the  analysis  will  be  given  here,  inasmuch  as  the  conventional 
staggered  procedure  is  not  recommended  for  general  use.  Only  one  illustrative  example  and 
some  general  observations  will  be  offered. 

The  stability  region  of  the  PE  formulation  is  always  x-6ounded,  i.e.,  parameters  u 
and  p play  no  role  in  the  analysis  (in  more  precise  terms,  the  worst  combination  is 
w = p = 0)  , and  neither  does  the  fluid  integrator  (8b).  For  a specific  structural  inte- 
grator (8a)  and  the  one-parameter  predictor  (13),  the  stability  regions  can  be  conveniently 
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displayed  in  the  (y>x)  plane.  Figure  1 shows  those  curves  for  the  trapezoidal  rule.  The 
largest  stable  x is  4,  which  can  be  attained  for  the  one-pass  solution  (k  = 1)  if 
Y = (the  predicted  pressure  value  is  the  mean  of  the  last  two  values).  If  the  process  is 

iterated,  the  peak  of  the  stability  region  is  reduced;  and  as  k " the  stability  limit 
approaches  the  predictor-independent  value  (18). 

It  should  be  stressed  that  a time  increment  limited  by  a x = h/5  of  order  2 to  4 is 
unacceptable  for  a general-purpose  integration  package.  The  constraint  can  be  practically  met 
only  in  a limited  class  of  problems,  such  as  ..ne  early-time  shock  response  of  submerged  struc- 
tures. If  the  response  is  dominated  by  low-frequency  structural  motions,  however,  that  incre- 
ment limit  is  intolerably  small;  so  small,  in  fact,  that  computational  error  accumulation  and 
computation  time  become  critical.  (The  largest  stable  x found  by  the  authors  is  8,  which 
requires  the  implementation  of  a fairly  elaborate  time-advancement  procedure.) 

3.3  OBSERVATIONS 


The  detailed  study  of  the  range  of  convergence  of  the  iterated  PE  formulation  can  be 
justified  by  the  following  considerations.  Suppose  we  had  found  that  the  iteration  (11)  con- 
verges for  the  entire  range  (or  at  least  a wide  range)  of  parameters.  The  converged  solution 
x(°°)  , ^(")  wouu  then  be  identical  to  that  furnished  by  the  f ully-implicit  scheme  resulting 


Figure  1 Stability  Limit  of  the  PE  Formulation  (11)  Treated  by  the 
Trapezoidal  Rule  and  the  Two-Step  Predictor  (13)  as  a 
Function  of  the  Number  of  Passes  (k)  per  Time  Step 
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Figure  2 Effect  of  Iterations  on  Stability 

Characteristics  of  Staggered  Procedures 

from  the  simultaneous  solution  approach,  as  sketched  in  Figure  2.  Under  such  assuptions, 
the  stability  region  of  the  iterated  PE  formulation  must  approach  that  of  the  f ully-implicit 
scheme.  The  latter  can  be  made  unconditionally  stable  by  selecting  a suitable  A-stable  LMS 
method.  Iteration  would  then  provide  a simple  strategy  for  transmuting  the  one-pass,  con- 
ditionally stable  semi-implicit  procedure  into  an  unconditionally  stable  scheme  of  greater 
accuracy;  this  strategy  might  have  led  to  a potentially  favorable  tradeoff  between  iteration 
cost  and  the  ability  to  utilize  larger  time  steps. 

As  increasing  k does  not  actually  improve  stability  to  any  significant  extent,  it 
appears  that  there  is  no  point  in  iterating  at  all  if  the  PE  formulation  is  used.  (In 
linear  problems,  an  iteration  cycle  costs  essentially  the  same  as  one  advancing  step,  but 
higher  accuracy  can  be  generally  achieved  through  equivalent  reductions  in  the  stepsize.) 

Finally,  it  should  be  noted  that  PE  iteration  can  be  made  to  converge  over  the  entire 
feasible  domain  of  parameters  by  the  application  of  a Jacobi  acceleration  strategy  [4 , § 4.3.1) 
to  (11).  The  resulting  iteration  process  is  virtually  identical,  however,  to  those  ensuing 
from  stabilization  of  the  original  equations  (4)  by  augmentation,  as  described  in  later 
sections . 
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3.4  VELOCITY  EXTRAPOLATION  PROCEDURE 

An  alternative  formulation  of  the  conventional  staggered  solution  procedure  can  be 
based  on  velocity  extrapolation  (VE)  applied  to  the  fluid  equation  (9b).  The  algorithmic 
properties  of  this  formulation  are  identical  to  those  of  the  PE  formulation,  and  need  not 
be  discussed  further. 


SECTION  IV 


f 


SOURCE  OF  INSTABILITY 

This  section  probes  more  deeply  into  the  underlying  causes  of  the  iteration  divergence 
and  temporal  stability  limitations  exhibited  by  the  conventional  (PE,  VE)  formulation  of 
the  staggered  solution  procedure.  First  we  examine  the  root  locus  diagram  of  the  character- 
istic equation  pertaining  to  the  f ully-implicit  formulation  of  the  model  system  (4),  noting 
its  essentially  stable  character.  The  PE  solution  procedure  is  then  cast  into  a 
differential-difference  (DD)  system  that  abstracts  the  time  integration  method.  The  analy- 
sis of  the  characteristic  equation  of  the  DD  system  shows  clearly  that  the  x-liroited 
instability  is,  as  expected,  caused  by  the  feedback  effect  of  the  extrapolated  pressure  term 
into  the  structural  equations  of  motion.  This  analysis  technique  has  the  important  advantage 
of  providing  results  that  are  independent  of  the  intrusion  of  a particular  time  integration 
scheme;  the  effect  of  the  latter  is,  in  fact,  of  secondary  importance.  Interpretation  of  the 
results  in  terms  of  a control  process  provides  direct  insight  into  techniques  for  the  stabil- 
ization of  the  staggered  solution  procedure  at  the  differential  equation  level;  such  tech- 
niques are  exploited  in  the  following  section. 

4.1  CHARACTERISTICS  OF  THE  FULLY- IMPLICIT  FORMULATION 

We  obtain  the  homogeneous  form  of  the  fully-implicit  model  system  by  transferring  the 
coupling  terms  in  (4)  to  the  left-hand  sides,  which  yields 

£ x + y + m2  x = 0 

(19) 

y-x+py  = 0 


Laplace  transformation  of  these  equations  without  consideration  of  initial  conditions  then 
yields 


_ 2 2 
C s + a 


s + p 


X(s) 

Y(s) 

(20) 


in  which  X(s)  and  Y(s)  denote  the  transforms  of  x(t)  and  y(t)  , respectively.  The 
associated  characteristic  equation  is  therefore 


(5  s2  + a2)  (s  + p)  + s2  = 0 


(21) 
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Figure  3 Root  Locus  Diagram  of  the  Fully-Xmplicit  System  (19) 
in  the  Laplace-Transform-Variable  Space 

The  root  locus  diagram  for  (21)  is  shown  in  Figure  3.  Because  all  physically  relevant 
branches  are  on  the  left-hand  side  (Re(s)  s 0,  (19)  are  inherently  stable  equations.  Con- 
sequently, a stable  numerical  solution  is  guaranteed  if  the  time  integrator  applied  to  (19) 
is  A-stable. 

4.2  CHARACTERISTICS  OF  THE  PE  FORMULATION 

A differential-difference  (DD)  equation  [5]  for  the  PE  formulation  can  be  obtained 
by  expressing  the  two-step  predictor  (13)  in  a delayed  continuous  form, 

yP(t)  = (1  + y)  y(t  - h)  - y y (t  - 2h)  (22) 


This  form  is  applied  to  the  model  system  (4)  to  obtain 


1 


which  is  the  DD  form  associated  with  the  single-pass  PE  formulation.  Note  that  (23)  is 
independent  of  any  time  integration  scheme.  Its  characteristic  equation  in  the  Laplace- 
transform  variable  s is  obtained  following  a procedure  identical  to  that  used  to  obtain  (21): 

a s2  + J)  (s  + M)  + S2  [(1  + y)  e~sh  - Y e~2sh]  = 0 (24) 


In  control  theory,  the  identifier  dead- time  element  is  often  used  for  terms  such  as  e 
— 2 s h 

and  e (see,  e.g.,  [6],  p.  284),  which  result  from  Laplace-transormation  of  time-delayed 

terms.  It  is  generally  acknowledged  that  the  occurrence  of  such  elements  in  control  loops 
has  a destabilizing  effect. 

In  the  limit  of  vanishing  stepsize  h,  (24)  approaches  the  fully  implicit  equation  (21), 
which  is  stable.  For  a finite  h,  it  can  be  shown  that  the  critical  combination  of  param- 
eters pertaining  to  the  stability  of  (24)  is  u>  « p = 0,  in  which  case  (24)  reduces  to 

C s + (1  + y)  e'sh  - Y e"2sh  = 0 (25) 

Equation  (25)  shows  that  stability  is  characterized  in  terms  of  three  parameters:  the 
"buoyancy  ratio"  £ , the  dead-time  constant  h (integration  stepsize),  and  the  extrap- 
olation parameter  y . The  critical  component  is,  of  course,  h , which  can  be  viewed,  from 
the  DD-form  standpoint,  as  the  sampling  rate  with  which  pressures  are  evaluated  and  fed  back 
into  the  structural  analyzer  by  the  extrapolator  (22).  The  destabilizing  effect  of  the  dead- 
time element  can  be  investigated  by  recasting  (25)  in  the  form 

1 + G(s)  = 0 (26) 


1 r x -sh  -2sh. 

— ( (1  + y)  e -Ye  ] 


The  most  expedient  approach  for  investigating  the  existence  of  unstable  roots 

Re(s)  >0  of  a complex  transcendental  equation  such  as  (26)  is  due  to  Nyquist  [7].  For 

the  special  form  (26),  Nyquist's  criterion  can  be  stated  as  follows:  let  the  Laplace- 

transform  variable  s encircle  the  entire  right-hand  plane  Re(s)  > 0 by  varying  s = 

2 

j v (j  = - 1)  from  v = - » to  + “ , and  find  the  angle  of  G(jv)  in  the  polar  co- 
ordinate system  [ j G (s ) | , £ G(s)]  that  satisfies 

* G(s)  = n , | G(s)  Is]  (2 


' 


Application  of  this  criterion  to  (27)  provides  the  two  conditions 


(1  + y)  cos  vh  - > cos  2vh  = 0 
E vh/h  z (1  + y)  si»  vh  - y sin  2vh 


(29) 


For  the  special  case  y = 0 , i.e.,  using  the  last  pressure  solution  as  the  predicted  value, 
(29)  gives  the  stability  limit 


X = h/(  s n/2  (30) 

Figure  4 shows  the  stability  limit  in  the  (x,  y)  plane  as  calculated  from  (29).  The 
destabilizing  effect  of  the  dead-time  element  is  now  apparent.  It  must  be  stressed  again 
that  the  stability  boundary  shown  in  Figure  4 applies  only  if  an  integration  scheme  of 
infinite  accuracy  (such  as  an  analytical  integration  or  a nontruncated  Taylor  series  algorithm) 
is  used  to  solve  the  system  (23),  and  that  the  boundary  is  perturbed  by  the  intrusion  of  a 
difference  time  integrator.  For  example,  the  use  of  the  trapezoidal  rule  as  the  structural 
integrator  increases  the  stability  domain,  as  can  be  deduced  from  a comparison  of  Figures  1 


INTERPOLATION  — - EXTRAPOLATION 


Figure  4 Stability  Limit  of  the  Difference-Differential  Form  (24) 
of  the  PE  Formulation  with  the  Two-Step  Predictor  (22) 
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and  A;  on  the  other  hand,  the  use  of  the  backward  Euler  method  (not  shown  here)  has  a detri- 
mental effect  on  stability. 

4.3  SUMMARY  OF  MAIN  FINDINGS 

• The  x-t>0UIKied  stability  of  the  conventional  (PE,  VE)  procedure  is  due  to 
the  dead-time  (delayed)  pressure  feedback  (22)  into  the  structural  equations 
(23).  This  process  can  be  viewed  (in  control  theory  terms)  as  a sampled- 
feedback  system,  in  which  the  pressure  feedback  is  delayed  by  a "hold" 

(dead-time  interval)  equivalent  to  the  stepsize  h.  This  representation  is 
illustrated  in  the  block  diagram  of  Figure  5. 

• The  most  commonly  used  strategy  for  stabilization  of  a system  that  includes 
dead-time  components  consists  of  introducing  a series  of  compensating 
elements  that  reduce  the  bandwidth  of  the  system's  frequency  response 
(see,  e.g.,  [8]  p.  118).  These  stabilization  techniques  essentially  amount 
to  the  introduction  of  damping  into  the  feedback  loop. 

• In  the  fully-implicit  solution  procedure,  damping  is  inherently  present  by 
virtue  of  the  energy  radiation  term  (y)  in  the  model  system  (19)  on  the 
left-hand  side.  The  effect  of  the  radiation  damping  on  the  structural 
response  roots  can  be  readily  appreciated  by  examining  Figure  3.  The 
structural  equation  (23)  of  the  PE  procedure,  however,  contains  no 
homogeneous  damping  term.  The  addition  of  artificial  damping  to  the 
structural  equation  (23a)  and/or  the  pressure  equation  (23b)  in  amounts 
sufficient  to  stabilize  the  solution  procedure  is  likely  to  have  cata- 
strophic effects  on  solution  accuracy.  We  are  therefore  motivated  to 
attempt  the  restoration  of  sufficient  damping  to  achieve  overall  stability 
through  appropriate  rearrangement  of  the  original  governing  equations. 


INCIDENT 

VELOCITY 


DEAD-TIME  ELEMENT 


Figure  5 Control  Loop  Representation  of  the  Pressure 
Extrapolation  Procedure  (24) 
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SECTION  V 


STABILIZATION 

The  foregoing  results  clearly  suggest  that  the  staggered  solution  procedure  may  be 
stabilized  by  the  addition  of  damping  terms  into  the  left-hand  side  of  the  governing 
equations  of  motion.  As  the  use  of  artificial  damping  is  ruled  out  by  accuracy  considerations, 
it  follows  that  the  governing  equations  must  be  tailored  in  such  a way  that  homogeneous 
damping  terms  either  appear  in  the  structural  equations,  or  be  added  to  the  fluid  equations, 
or  both.  We  now  proceed  to  describe  three  stabilized  formulations  generated  through  such 
"damping  augmentation"  techniques. 

9 

5.1  PRESSURE- INTEGRAL  EXTRAPOLATION  (PIE)  FORMULATION 

The  homogeneous  portion  of  the  structural  equations  for  the  f ully-implicit  solution 
procedure  (19)  effectively  exhibits  damping  due  to  the  simultaneous  velocity  (x)  feed- 
back from  the  coupled  pressure  equation  (19).  On  the  other  hand,  the  structural  equation 
for  the  PE  formulation  (23)  is  excited  by  a time-lagged  velocity  feedback;  a destabi- 
lizing effect  occurs  because  of  the  delayed  energy  dissipation  in  the  feedback  loop  (cf. 

Figure  5).  We  can  correct  this  situation  by  eliminating  the  time-lagged  velocity  feed- 
back; this  is  accomplished  mathematically  by  combining  (4a)  and  (4b): 


£ x + 


2 

u x 


y = - (x  - y y) 


(31) 


which,  upon  transferring  the  damping  terra  (x)  to  the  left-hand  side,  becomes 


£ X + X +10 


2 

x 


v y 


(32) 


The  solution  procedure  based  upon  (32)  along  with  the  original  pressure  equation  (4b)  will 
be  called  the  pressure- integral  extrapolation  (PIE)  formulation,  inasmuch  as  the  pressure- 
integral  y is  involved  in  the  prediction  process. 

5.2  DISPLACEMENT  EXTRAPOLATION  (DE)  FORMULATION 

Stabilization  can  also  be  achieved  by  augmentation  of  the  fluid  equation  so  as  to 
increase  the  pressure  decay  rate  in  the  homogeneous  system.  Qualitatively  speaking,  this 
strategy  works  as  long  as  the  stability  "margin"  of  the  modified  pressure  equation  over- 
whelms the  destabilizing  effect  of  the  structural  equations.  The  modified  pressure 
equation  is  obtained  by  the  substitution  of  (4a)  into  the  time-differentiated  (4b): 

y + y y = x = - (y  + u>2  x)/£  (33) 
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which  can  be  rearranged  to  form 


y + (u  + I/O  y 


x/5 


(34) 


The  solution  procedure  based  upon  (34)  along  with  the  original  structural  equation  (4a)  will 
be  called  the  displacement  extrapolation  (DE)  formulation. 

5.3  OTHER  FORMULATIONS 

The  three  formulations  of  the  model  problem  described  so  far  are  collected  in  Table  2. 
It  turns  out  that  additional,  more  complicated,  formulations  may  be  constructed.  As  these 
have  not  been  found  to  possess  any  particular  advantage  over  the  two  preceding  formulations, 
they  will  not  be  discussed  here. 


TABLE  2.  Model  Formulations  for  Staggered  Solution  Procedures 


SECTION  VI 


ITERATION  CONVERGENCE  OF  STABILIZED  FORMULATIONS 

The  iteration  convergence  properties  of  the  stabilized  formulations  derived  in  the 
foregoing  section  may  be  studied  with  a technique  similar  to  that  used  to  analyze  the 
iterated  PE  scheme  (11).  The  convergence  conditions  are  collected  in  Table  3. 

It  is  easily  shown  that  the  spectral  radii  for  the  two  augmented  schemes  are  less  than 
unity  over  the  entire  feasible  domain  (0  to  + 00 ) of  the  parameters  x « ft  and  V . The  key 
significance  of  this  result  is:  if  the  solution  is  iterated  to  convergence  at  each  time 
step,  the  solution  of  the  fully- implicit  scheme  is  recovered  (cf.  Figure  2).  With  such  con- 
vergence guaranteed,  we  have  now  at  our  disposal  a feasible  computational  strategy  by  which 
unconditional  stability  can  be  attained  within  the  basic  organization  of  the  staggered  solu- 
tion procedure. 


TABLE  3.  Iteration  Convergence  Characteristics  of  Staggered  Solution  Procedures 


Formulation 

Spectral  Radius,  k , of  Iteration  Matrix 
Iteration  converges  if  k < 1. 

PE 

ex  x/Ui  + 82  ft2)  (1  + By  ¥)] 

PIE 

By  Bx  x */[(l  + 8y  ¥)  (1  + Bx  X + 82  ft2)] 

DE 

6 B2  x ft2  /Id  + 8 ¥ + 6 x)  ( 1 + 62  ft2)] 

y x y y x 

6.1  SELECTION  OF  A PREDICTOR 

A wide  variety  of  predictor  formulas  can  be  used  in  the  practical  implementation  of  the 
stabilized  procedures.  The  ideal  formula  should  be  able  to  provide  satisfactory  extrapola- 
tion accuracy  for  physically  meaningful  response  components  without  jeopardizing  the  temporal 
stability  characteristics;  this  is  indeed  a delicate  compromise.  The  following  considera- 
tions are  aimed  at  restricting  the  class  of  admissible  formulas  to  predictors  that  exclude 
"historical  derivatives".  A predictor  that  includes  derivatives  of  past  solutions  intro- 
duces large  extrapolation  errors  in  both  high-frequency  structural  response  components 

(ft  > > 1)  and  rapidly-decaying  pressure  response  components  (¥  > > 1)  . Table  3 shows 

2 

that,  as  either  x ^ or  X ft  becomes  large,  the  corresponding  spectral  radii  approach 
unity.  On  the  other  hand,  for  low-frequency  components,  i.e.,  as  either  x or  ¥ and/or 
ft  tend  to  zero,  all  spectral  radii  approach  zero.  This  means  that  extrapolation  errors 
committed  for  "noise"  components  of  the  response  (i.e.,  those  with  large  ¥ and/or  ft  ) will 
decrease  very  slowly  with  increasing  iteration  index. 
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We  therefore  restrict  our  considerations  to  predictors  that  are  less  prone  to  excite 
noise  components.  The  following  three-step,  two-parameter  family  can  be  considered  suf- 
ficient for  our  envisioned  applications: 

Zn  = ((1  + V (1  ' Y2}  + 3 V Vl 

- Id  - Y2)  Yj  + 3 y2]  zn_2  + y2  zn_3 

where  z stands  for  the  solution-state  term  being  extrapolated.  Special  cases  of  (35)  are 
listed  in  Table  4. 

6.2  SINGLE-PASS  STABILIZED  PROCEDURES 

The  difference  equations  for  the  single-pass  PIE  formulation  are  [cf.  (T2 ) ] 

E x(1)  = B2  x f y(0)  + (1  + 0 x)  hx  + 6 h hx 
xn  x n x n x n 


E y(1)  = 0 (x(1>  - hX)/0  + hy 

y n y n n x n 


where  y^  ' is  predicted  by  means  of  (35), 


2 2 

E =1  + 3 X + B n 
x x x 


TABLE  4.  Predictors  Considered  for  the  Present  Study 


Case. 

Parameters 

Order  of  Accuracy 

I 

= Y2  = 0 

Zero  order 

IX 

= 1/2,  Y2  - 0 

Improved  zero-order 

III 

Yi  = i.  y2  = 0 

First  order 

IV 

Y i = 1.  Y2  = 1/2 

Improved  first  order3 

V 

Y’l  = 1,  Y2  = -2/3 

Lease-square  fit 

aThis  formula  is  the  basis  of  the  widely  used  Houbolt  integrator  (9) 

1 


and  Ey  is  given  in  (12b) . The  temporal  stability  of  (36)  may  be  examined  by  seeking  a 
nontrivial  solution  of  the  form 


where 


z = X z . 

— n — n-1 


z = (x,  x,  x.  i) 


and  requiring  a bounded  solution  for  (36),  viz. 


X s 1 


(38) 


(39) 


Difference  equations  similar  to  (36)  may  also  be  obtained  for  the  DE  formulation.  Evaluation 
of  the  stability  condition  (39)  for  various  time-integrator/predictor  combinations  leads  to 
tedious  algebraic  manipulations  (Appendix  C) . Table  5 summarizes  the  results  for  the  single- 
pass implementation  of  the  formulations  listed  in  Table  2. 


TABLE  5.  Stability  Characteristics  of  Single-Pass  Solutions 


Formulation 

Stability  Properties 

Pressure  Extrapolation  (PE) 

X < 4 for  trapezoidal  rule3 

Pressure  Integral 

Unconditionally  stable*3  for 

Extrapolation  (PIE) 

predictor  cases  I,  IIC 

Displacement 

Unconditionally  stable*3  for 

Extrapolation  (DE) 

predictor  cases  I,  III  when 
trapezoidal  rule  is  used^ 

aSee  Figure  1. 

^"Unconditionally  stable"  means  temporal  stability  for  any  feasible 
combination  of  physical  parameters  and  integration  stepsize. 

CSee  Table  4. 

^Allowable  predictors  depend  on  integration  method  used. 
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6.3  OBSERVATIONS  ON  SINGLE-PASS  SOLUTION  STABILITY 


The  single-pass  stability  properties  listed  in  Table  5 account  for  the  combined  effects 
of  the  predictor  and  the  time-integrator.  The  results  indicate  that  the  two  augmented  formu- 
lations are  capable  of  delivering  stable  single-pass  solutions  provided  that  the  extrapola- 
tion accuracy  is  suitably  restricted.  Such  restrictions  are,  in  general,  affected  by  the 
integration  method.  For  example,  application  of  the  trapezoidal  rule  to  the  DE  formulation 
limits  the  extrapolator  accuracy  to  first  order  (case  IV  of  Table  4),  whereas  the  application 
of  an  A-stable  backward-difference  scheme  allows  the  use  of  all  of  the  predictors  listed  in 
Table  4. 

We  recall  that  the  two  stabilized  formulations  can,  if  iterated  to  convergence,  attain 
the  unconditional  stability  of  the  fully-implicit  solution  procedure  (cf.  Table  3 and  Figure 
2),  regardless  of  the  choice  of  extrapolator.  This  observation  raises  again  the  issue  of 
finding  a cost-minimization  compromise  between  a single-pass  solution  process  with  small 
time  increments  and  an  iterated  solution  process  with  larger  time  increments.  A resolution 
of  this  issue  must  await  the  accumulation  of  experience  in  the  application  of  these  new 
techniques  to  actual  linear  and  nonlinear  fluid-structure  interaction  problems. 

Finally,  it  should  be  mentioned  that  stability  analyses  of  the  differential- 
difference  equations  pertaining  to  the  single-pass  PIE  and  DE  formulations  were  also  per- 
formed using  integral-transform  techniques  similar  to  those  used  for  the  PE  formulation. 

These  analyses  provided  valuable  insight  into  the  sensitivity  of  the  stability  character- 
istics to  the  various  parameters  appearing  in  the  model  formulations.  The  results  will  not 
be  reported  here,  however,  because  of  the  length  and  complexity  of  the  algebraic  manipula- 
tions involved. 


SECTION  VII 


NUMERICAL  EXPERIMENTS 

A series  of  numerical  experiments  was  conducted  with  dual  objectives:  to  verify  the 
predictions  of  the  convergence  and  stability  analyses,  and  to  assess  the  global  accuracy 
of  the  staggered  solution  procedures.  We  have  adopted  a direct  accuracy  assessment  over  an 
analytical  accuracy  analysis,  as  the  former  provides  an  overall  (global)  measure  of  accuracy, 
which  results  from  the  combined  effects  of  integration  algorithm  and  implementation  form  [10], 
as  well  as  extrapolation  process.  The  accuracy  measure  chosen  for  the  present  study  is 


e 


i=l 


1/2 


(40) 


where  the  superscripts  ( ) and  ( ) denote  the  exact  and  computed  solutions  at  the  dis- 
crete time  station  t = t^  , respectively,  and  N is  the  total  number  of  time  steps  taken 
to  evaluate  e . The  exact  solutions  were  obtained  by  using  the  characteristic  roots  of 
the  f ully-implicit  procedure  (21)  and  applying  the  initial  conditions 


x(0)  = (£  x/Q)2 

x(0)  = 1 

y(0)  = £ 


These  values  produce  reasonably  scaled  response  histories  for  all  three  response  variables 
of  interest. 

A number  of  implicit,  one-derivative,  A-stable  LMS  methods  were  used  for  the  numerical 
solution  of  each  of  the  (4),  among  them  being  the  trapezoidal  rule  and  the  backward- 
difference  operators  of  Park  [11]  and  Gear  [12].  In  the  results  to  be  presented  here,  the 
trapezoidal  rule  was  used  to  advance  the  structural  response  while  the  Park  three-step 
method  was  applied  to  the  solution  of  the  fluid  equation. 

Error  comparisons  after  20  steps  were  made  for  the  structural  displacement  x(t)  and 
velocity  x(t)  as  well  as  the  fluid  pressure  y(t)  ; however,  only  the  displacement  results 
will  be  shown.  Accuracy  assessments  were  carried  out  for  the  parameter  range 

1 s X « 10,  0-01  « 0 s 1.0,  0.0001  s T s 0.1, 

with  the  results  being  relatively  insensitive  to  V . For  this  reason,  comparisons  are 
shown  here  only  for  T = 0.1,  which  value  usually  produced  the  greatest  error. 
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Figure  6 shows  a comparison  between  the  pressure  integral  extrapolation  scheme,  denoted 

PIE  DE 

by  x , versus  the  displacement  extrapolation  scheme,  denoted  by  x .In  each  case  a 

two-step  standard  linear  extrapolation  method  (case  III  of  Table  4)  was  used.  Although  the 

PIE  solution  displays  superior  accuracy,  this  scheme  violates  the  practical  requirement 

(justified  in  the  following  section)  that  the  structural  equations  be  left  unchanged.  It  is 

therefore  important  to  show  that  the  accuracy  of  the  DE  scheme  can  be  improved  without 

incurring  a substantial  increase  in  computation  time.  Two  immediate  ways  of  accomplishing 

this  objective  are:  the  use  of  a three-step  extrapolation  method  (case  IV  of  Table  4), 

and  the  performance  of  one  iteration  on  the  fluid  equation  solution.  The  latter  effectively 

amounts  to  an  additional  half-pass  through  the  solution  strategy. 

As  the  integration  of  the  structural  equations  normally  controls  the  computation  time 
in  large-scale  problems,  this  additional  iteration  can  be  carried  out  with  a modest  increase 
in  the  total  run  cost. 

The  resulting  improvement  in  the  accuracy  of  the  DE  procedure  is  shown  in  Figures  7 and 
8,  where  the  first  subscript  indicates  whether  a two-step  or  a three-step  extrapolation 
method  is  used  and  the  second  subscript  denotes  the  number  of  iterations  made  on  the  fluid 
equation.  Hence  xl?E  is  identical  to  xDE  in  Figure  6.  As  can  be  seen,  the  solution 

4 , U 

accuracy  is  more  sensitive  to  the  extrapolation  scheme  used;  however,  the  value  of  applying 
one  iteration  to  the  fluid  equation  is  also  apparent.  Similar  accuracy  improvements  were 
observed  for  the  velocity  and  pressure  components  of  the  solution. 


STRUCTURAL  CRI'XrlNlr  0 


Figure  6 Global  Error  for  PIE  and  DE  Forms 
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SECTION  VIII 


IMPLEMENTATION  CONSIDERATIONS 

Thus  far  we  have  concentrated  on  stabilization  of  the  staggered  solution  procedure;  we 
have  succeeded  in  developing  several  stabilized  forms  of  the  two  degree-of-f reedom  system  (4). 
Application  of  the  two  formulations  of  greatest  interest  to  the  original  matrix  system  (1) 
provides  the  equations  shown  in  Table  6.  Before  selecting  one  of  these  formulations  for 
implementation  in  a fluid-structure  analysis  code,  it  is  necessary  that  we  discuss  the 
practical  implications  associated  with  such  a decision. 

If  there  were  no  a priori  constraints  with  regard  to  software  development,  i.e.,  the 
program  developer  had  complete  freedom  to  construct  both  the  structural  and  fluid  analyzers 
from  scratch,  the  overriding  selection  factor  would  probably  be  the  accuracy  characteristics 
displayed  by  each  formulation.  As  it  is,  there  now  exist  many  large-scale  linear  and  non- 
linear structural  analyzers  that  incorporate  capabilities  for  transient  response  analysis. 


TABLE  6.  Matrix  Implementation  Forms  for  Stabilized  Procedures 


Formulation 

Implementation 

Pressure  Integral 
Extrapolation  (PIE) 

MU  + (D  + p cTATT)u  + Ku=f^  + r 

- T A (s1  - p c u1)  +p  c T AM^1  A a 

. -1  T I 

A q + P c A Mf  A £ = p c A (T  u-u) 

Displacement 
Extrapolation  (DE) 

Mii+Du+Ku  = f +r-TA(q  + pI) 

A £ + p c (A  Hj'  A + A M G A)  £ = - p cAu* 

. m-G  . I 
- p c A M A £ 

+ p c A M_G  C1TT(f  +r-Du-Ku) 

where  M_G  = CT  (TT  M T)"1  C , C = TT  T 

Inspection  of  Table  6 reveals  that  only  the  DE  formulation  has  no  impact  upon  the 

structural  equations  of  motion,  while  the  PIE  formulation  (as  well  as  the  other,  more  com- 

T 

plicated,  formulations)  requires  the  addition  of  a damping  term,  p c T A T , to  the  left- 
hand  side  of  the  structural  equation.  This  term  can  be  expected  to  have  an  adverse  effect 
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upon  the  sparseness  characteristics  of  the  dynamic  coefficient  matrix  [the  multi-degree- 

of-freedom  analog  of  E in  Eqs.  (11)  and  (36)],  because  each  fluid  boundary  element  on 

X t 

the  contact  or  "wet"  surface  B will  normally  overlap  several  structural  elements  . Addi- 
tional nonzero  coefficients  then  appear  outside  the  "profile"  of  the  structural  stiffness 
matrix  in  entry  positions  pertaining  to  the  wet  structural  displacements  on  B.  As  matrix 
connectivity  characteristics  are  generally  controlled  by  the  structural  grid  information, 
extensive  and  expensive  software  modifications  would  be  required  to  include  such  a damping 
tern  in  existing  "dry-structure"  analyzers. 

For  the  analysis  of  linear  structures,  a general-purpose,  stand-alone  integration 
package  based  on  available  A-stable  LMS  methods  could  be  built  for  the  PIE  and  other 
formulations,  which  would  accept,  as  input  data,  preprocessed  matrices  assembled  by  existing 
linear  structural  analysis  codes  (with,  perhaps,  some  restrictions  on  sparse-matrix  storage 
formats).  However,  for  structures  exhibiting  nonlinear  behavior,  such  an  undertaking  would 
necessitate  the  coupling  of  the  nonlinear  solver  with  the  dedicated  integration  package  to 
provide  continuously  updated  information  in  the  form  of  factored  matrices,  pseudo-force 
vectors,  and  the  like.  In  this  case,  the  driving  consideration  would  certainly  be  preser- 
vation of  the  autonomy  of  the  structural  analyzer,  a restriction  that  precludes  the  use 
of  either  the  PIE  or  other  forms. 

An  additional  argument  that  reinforces  this  conclusion  is  the  fact  that  the  DAA  is 
only  the  lowest-order  member  of  a family  of  surface  interaction  approximations  [3].  The 
next  member  of  that  family  is  characterized  by  a second-order  matrix  ODE  that  would  replace 
(lb).  The  inclusion  of  a higher  derivative  of  the  scattered  pressure  has  the  effect  of 
adding  widely-connected  matrix  terms  to  both  the  damping  and  stiffness  matrices  in  the 
structural  equations  if  the  equivalents  of  the  PIE  and  other  formulations  are  adopted. 

In  summary,  it  seems  clear  that  the  displacement  extrapolation  formulation  is  the 
only  stabilized  procedure  that  satisfies  the  practical  requirements  of  modularity  and 
avoids  the  proliferation  of  special-purpose  versions  of  existing  structural  analyzers. 


The  fluid  discretization  can  be  much  coarser  than  the  structural  discretization  because 
the  calculation  of  spatial  displacement  gradients  is  only  required  for  the  structural 
model  in  order  to  obtain  strains  and  stresses. 
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SECTION  IX 


CONCLUDING  REMARKS 

The  original  goal  of  this  investigation  was  to  find  an  optimal  implementation  of  the 
staggered  solution  strategy  for  the  structure/DAA  equations,  with  the  model  equations (3) 
or  (4)  as  a point  of  departure.  It  was  hoped  that  such  an  implementation  would  display 
adequate  stability  and  accuracy  properties  to  provide  a suitable  alternative  to  the  fully- 
implicit  (simultaneous-solution)  approach,  but  without  the  computational  drawbacks  of  the 
latter.  Although  that  goal  was  attained  beyond  our  expectations,  the  study  guidelines  had 
to  be  frequently  redefined  along  the  way. 

It  is  a common  attribute  of  intricate  research  projects  that  essential  ingredients  of 
the  problem,  i.e.,  those  aspects  crucial  to  success  or  failure,  are  seldom  recognized  in 
advance.  The  work  reported  here  is  no  exception.  An  extensive  initial  series  of  parametric 
studies  aimed  at  extending  the  stability  limit  of  the  conventional  staggered  solution 
procedure  failed  to  yield  satisfactory  results.  When  attempts  to  introduce  artificial 
damping  terms  proved  fruitless,  and  when  we  were  seriously  considering  the  adoption  of 
explicit  integration  methods,  the  augmentation  concept  was  tried,  which  led  to  the  deri- 
vation of  the  PIE  formulation.  That  serendipitous  discovery  prompted  a systematic  develop- 
ment of  other  formulations  and  the  final  selection  of  the  displacement  extrapolation 
formulation. 

The  interpretation  of  staggered  solution  procedures  from  the  viewpoint  of  control 
theory  emerged  in  final  form  as  this  report  was  being  prepared.  This  interpretation  is 
deemed  important  for  a variety  of  reasons: 

• The  effects  of  the  predictor  can  be  incorporated  into  differential-difference 
equations,  thus  isolating  the  effects  of  the  integration  method.  (The  latter 
can  in  fact  be  included,  if  desired,  by  adjoining  to  the  characteristic 
determinant  terms  resulting  from  a Laplace-transformation  of  the  integration 
operator  [13].) 

■ A large  body  of  control  theory  techniques,  such  as  frequency  bandwidth  reduction, 
damping  augmentation,  etc. , can  be  applied  to  the  stabilization  of  the  system  in 
transform  space,  and  such  remedies  reverted  to  the  time  domain  and  interpreted 
in  terms  of  computational  procedures;  see,  for  example,  [14,  15]. 

• A block  representation,  such  as  the  one  shown  in  Figure  5,  provides  a convenient 
framework  for  setting  up  parameter  studies  using  a digital  or  analog  computer. 

• The  interpretation  is  applicable  to  general  classes  of  coupled  field  problems, 
rather  than  being  confined  to  fluid-structure  interaction  analysis. 
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With  regard  to  the  last  of  these,  there  is  presently  a growing  interest  in  the  application  of 
computerized  analysis  methods  to  a host  of  problems  that  are  modeled  through  coupled-field 
evolutionary  equations  arising  in  geomechanics,  biomechanics,  thermoelasticity  and 
magnetohydrodynamics,  to  cite  only  a few.  Staggered  solution  procedures  appear  parti- 
cularly attractive  when  software  modules  (analyzers)  are  available  for  processing  the 
individual  uncoupled  problems.  These  modules  can  be  connected  to  form  a serially- 
executable  "analyzer  network"  through  interfacing  mechanisms  based  on  the  prediction  of 
appropriate  subsets  of  the  complete  solution  vector.  Stabilization  of  these  mechanisms 
at  the  governing  ODE  level  is  most  effective  for  avoiding  "delayed  feedback"  instability. 

Once  satisfactory  stabilized  formulations  are  synthesized,  accuracy  and  implementation 
considerations  can  be  used  to  select  the  most  desirable  formulation  and  associated  extra- 
polation formulas. 
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APPENDIX  A 


1 f 

. 


DERIVATION  OF  MODEL  EQUATIONS 
Consider  the  following  form  of  the  homogeneous  system  (2) 


M ii  + K u 


= -T  A q 


(A-l) 


Aq  + pcBq  = 


T . 
CAT  u 


in  which 


B = A M, 
f 


-1 


A . Equations  (A-l)  are  obtained  upon  setting  D = 0 (no  structural 

T —1  —1 

damping)  and  premultiplying  the  fluid  equation  by  A Mf  = A ; the  latter  operation 
makes  the  coupling  matrices  in  the  right-hand  side  of  (A-l)  the  transpose  of  one  another. 
The  symmetric  algebraic  eigenproblems  associated  with  the  left-hand  side  of  (A-l)  are 


<-“i  % + V «i 


0,  i = 1,  . . . n 


(A-2) 


< 4 * B)  - 0,  J - 1. 


The  u±  are  the  dry  structural  modes  (also  called  in  vacuo  modes)  pertaining  to  the 
undamped  natural  frequencies  . The  q^  are  fluid  boundary  modes  associated  with  the 
pressure-decay  exponents  = p i . (These  fluid  modes  are  discrete  analogs  of  the 
eigenfunctions  of  the  added  mass  tensor  associated  with  the  Laplace  equation,  subject  to 
appropriate  orthonormality  conditions  with  respect  to  surface  area  [16].)  In  accordance 
with  (A-2),  we  introduce  the  transformations 


u = Sw;  q = F g 


(A-3) 


where  the  matrices  S and  F are  formed  with  the  dry-structure  and  fluid-boundary  modes, 
respectively,  lined  up  as  columns  of  unit  length,  and  the  vectors  w and  g collect  the 
associated  mode  amplitudes  (generalized  coordinates) . The  transformed  system  (A-l)  is  then 


m w + k w 
— ~s  — 


a g + pcjfcg 


~ S I A £ S ‘ -H  .8 


T T T ■ 

pcF  AT  Sw  = pcH  w 


(A-4) 


where  m^  , , a and  jj  are  diagonal  matrices  of  generalized  quantities: 


A-l 


m,  = S MS,  = S' 1 K S 


a = F A F 


(A-5) 


F B F = 


a <IT  £>  1 a = a Ef1  a 


Although  the  left  side  of  (A-4)  is  uncoupled,  each  dry  structural  mode  will,  in  general, 
be  coupled  to  each  fluid  boundary  mode  through  the  right-hand  side  terms.  Consequently,  for 
each  mode  pair  u^  , , we  obtain  a two  degree-of-f reedom  system  in  the  associated  gen- 


eralized coordinates  w^  and  g^  : 


m.w.+k.w.  = -<?.  . a.  g. 

si  1 si  l lj  J J 


a.  gj  + P c b.  g.  = pea.  9jl  wt 


(A-6) 


In  (A-6),  the  <p.  . are  modal  coupling  coefficients  between  the  i-th  structural  mode  and 
the  j-th  fluid  boundary  mode,  i.e.. 


hij  2 A Sj  Ei  ~ q-i 


'ij 


j T A T A 

q . A q . q . A q . 

-J  ~ -J  -J  ~ -J 


where 


(A-?) 


(A-8) 


r.  = u?  T . 

-i  -i  ~ 

We  will  now  show  that  | <?_  | ^ 1 if  all  entries  of  the  structure-to-f luid  grid 
transformation  matrix  are  nonnegative  (a  restriction  that  can  always  be  met  by  using 
appropriate  "area-lumping"  schemes).  Equilibrium  of  forces  normal  to  the  contact  surface 
demands  that  the  sum  of  entries,  °f  each  column  t^  of  T be  unity: 


£ tk4  ^ I Sa  ' ” 1 

k=l  k 


(A-9) 


Hence 


1,  »2  ' 


(A-10) 


Using  the  norm  inequality  II  x y II  z II  x II  • II  y II  , we  can  easily  show  that  II  II  2 1 • 


A- 2 


V 0 


Now  from  (A-7)  we  obtain 


(ij  % ~ Ei>  A q.  =0 


(A-ll) 


As  A is  positive  definite  and  c^.  is  not  identically  zero,  it  follows  that  the  vector  in 

parenthesis  must  be  orthogonal  to  A q . Premultiplying  that  vector  by  and  solving 

for  V.  . , we  obtain  the  desired  result 
rj 


<P.  • 

ij 


q . r . 
-J  -r 
T 

q-J  q-J 


ij  Ei  ^ 1 


(A-12) 


Eq.  (A-12)  indicates  that  setting  <P„  =1  in  (A-6)  implies  a conservative  assump- 
tion on  the  modal  coupling,  which  is  the  essential  source  of  stability  difficulties  in  the 
staggered  solution  procedures.  The  worst  modal  coupling  case,  i.e.,  <p  „ = 1 , occurs  when 
both  structural  and  fluid  boundary  modes  coincide,  as  in  the  case  of  breathing  motions  of  a 
submerged  spherical  shell.  Assuming  identical  surface  grids  for  this  case,  we  have 


and  = 1 
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SUBMERGED  SPHERICAL  SHELL 

To  get  an  idea  of  the  range  of  values  assumed  by  the  various  parameters  that  appear 
in  the  study  of  the  staggered  solution  procedures,  it  is  instructive  to  study  a closed,  thin 
spherical  shell  submerged  in  an  infinite  fluid.  The  dry-structure  and  fluid-boundary  modes 
for  the  continuum  problem  coalesce  and  can  be  analytically  expressed  in  terms  of  Legendre 
polynomials  in  cos  6 , 0 being  the  meridional  angle.  An  analysis  of  the  modal  spectrum 

[17  ] pertaining  to  axisymmetric  modal  motions  with  n circumferential  waves  gives  for  the 
generalized  structural  mass,  stiffness,  fluid  mass  and  contact  areas. 


m 

sn 


4irp  tr2  [1  + n (n  + 1)42  ] / (2  n + 1) 
s n 


k 

sn 


_2 

m u) 
sn  sn 


m = 4irp  r3  / [(n  + 1)  (2n  + 1)] 

a = 4ir  r2  / (2n  + 1) 

n 


(B-l) 


where  pg  , t and  r are  the  shell  density,  thickness  and  radius,  respectively, 

denotes  modal  natural  frequency  [not  to  be  confused  with  the  reduced  frequency  u)  defined  in 

(5b)]  and  5 are  ratios  of  tangential  to  radial  modal  amplitudes,  given  by  [18,  Ch.  10] 

4 = (1  + v)  / [-w2  + n2  + n + (1  - v)]  ( B— 2 ) 

n n 

where  v is  Poisson's  ratio,  and  the  are  dimensionless  in-vacuo  natural  frequencies 
normalized  to  cg/r  > cs  being  the  plate  velocity  for  the  shell  material).  The 
dimensionless  parameters  £ , w and  p can  be  obtained  from  Eqs.  (5),  in  which  l = r , 
as 


= [1  + n (n  + 1)  C2  ] pgt/pr) 


0)  = £ U)  C /C 

n ns 


p = n + 1 
n 

Values  of  ( , u and  p are  compiled  in  Table  7 for  a steel  shell  in  water  (cs/cs3.53, 
pg/p  = 7.67,  v = 0.30)  with  a thickness-to-radius  ratio  of  0.01.  The  lower  (upper) 
eigenfrequency  branch  embodies  predominantly  flexural  (membrane)  motions. 


B-l 


I 

I 


Note  that  the  lowest  values  of  these  parameters  are  determined  by  physical  and  geometric 
characteristics  of  the  fluid-structure  problem,  whereas  the  high  values  are  determined  by 
truncation  of  the  modal  spectrum,  lor  the  discrete  model,  the  largest  values  are  essentially 
functions  of  the  mesh  fineness.  These  considerations  should  be  helpful  in  the  order-of- 
magnitude  assessment  of  the  parameter  value  range  for  more  complex  problems. 


TABLE  7.  Values  of  £ , to  , and  p for  Some  Axisymmetric 
Modes  of  a Steel  Shell  in  Water  with  t/r  = 0.01 


Mode  Order 

Branch  n 

Cn 

A 

co 

n 

z 

n 

CO 

n 

Mn 

Lower  1 

1.000 

0 

0.230 

0 

2 

2 

0.270 

0.701 

0.110 

0.821 

3 

3 

0.123 

0.830 

0.091 

0.882 

4 

n»l 

-2 
~ n 

~l-l/2n 

-0.077 

-0.977 

n + 1 

Upper  0 

0 

1.61 

0.077 

0.977 

i 

1 

-0.500 

1.98 

0.115 

2.37 

2 

2 

-0.616 

2.72 

0.251 

4.81 

3 

3 

-0.680 

3. 64 

0.502 

9.  10 

4 

n»  1 

~-n/n+l 

~(n2  + 3)** 

~0.077n2 

~0.98n2 

n + 1 
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A STABILITY  ANALYSIS  OF  THE  PIE  FORMULATION 


Let  us  write  (6)  and  (7)  as 


p (Z  -k)  - 6o(Z  , ) = 0,  k = 1 , - - , m 

~n  ~n-k 


» «n-k>  'E'tVl  •««„-»>  "Z 


8k  Vk 


and  the  PIE  model  equations  in  Table  2 as 


0 10 


-10  0 


C00  Zn  + 1 ui  -u 


0 0 1 


0 0 p 


where  Z = [x,  x,  y] 

For  illustrative  purposes,  introduce  the  first-order  extrapolator  £see  (35)  and 
Table  4]. 


x ■ 2x  , - x _ 
n -n-1  ~n-2 


combining  (C-l)  through  (C-3) , we  obtain 


0 10 


0 0 1 


-10  0 


C 0 0 P(Zn_k)  + 6 1 ~v\  °'Zn-k)  = 5 0 


0 0 p 


0(in-k> 


Substitution  of  (38)  into  this  equation  then  yields 


0 1 0 


0 p(A)  + <S 


Z =0 
— n-m 


0 0 1 


r.O  p 


C-I 


where 


m 


m 

p(A)  = “k  Am  k » °(a>  = y]\Am~k 

k=0  k=0 

Now  a nontrivial  solution  of  (C-6)  exists  only  when 


1 

'0 

1 

o' 

'-I  0 

O' 

det  1 

c 

0 

0 

P (A)  + 5 

1 w2 

-y 

a(  A) 

0 

0 

1 

* 0 ’ 

y 

\ 

. 

. 

L A A 

from  which  we  obtain  the  following  characteristic  equation 


p3(X)  + S p p (A)  a (A)  rp (A)  - (A)  (2A 


D/A2] 


cf 3 ( A ) + -^Do2(A)  p(A)  + | p2(A)  o(A)  + ^ P(A>  °2 (A)  = 0 


(C-7) 


(C-8) 


(C-9)  _ 


Therefore, 
absence  of 
uh  , yh  , 
results  of 


the  search  for  unconditional  stability  is  reduced  to  determining  the  complete 
any  root  of  (C-9)  with  magnitude  greater  than  unity  for  all  possible  values  of 
and  h/C  in  combination  with  the  time  integrator  and  extrapolator  used.  The 
such  determinations  are  summarized  in  Table  5. 
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